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Abstract. We present the results of a study confronting density maps reconstructed by the Delaunay Tessellation 
Field Estimator (DTFE) and by regular SPH kernel-based techniques. The density maps are constructed from 
the outcome of an SPH particle hydrodynamics simulation of a multiphase interstellar medium. The comparison 
between the two methods clearly demonstrates the superior performance of the DTFE with respect to conventional 
SPH methods, in particular at locations where SPH appears to fail. Filamentary and sheetlike structures form 
I telling examples. The DTFE is a fully self-adaptive technique for reconstructing continuous density fields from 

, discrete particle distributions, and is based upon the corresponding Delaunay tessellation. Its principal asset is 

fT^ ' its complete independence of arbitrary smoothing functions and parameters specifying the properties of these. 

I As a result it manages to faithfully reproduce the anisotropics of the local particle distribution and through its 

CO . adaptive and local nature proves to be optimally suited for uncovering the full structural richness in the density 

' distribution. Through the improvement in local density estimates, calculations invoking the DTFE will yield a 

much better representation of physical processes which depend on density. This will be crucial in the case of 
^ feedback processes, which play a major role in galaxy and star formation issues. The presented results form an 

encouraging step towards the application and insertion of the DTFE in astrophysical hydrocodes. We describe an 
outhne for the construction of a particle hydrodynamics code in which the DTFE replaces kernel-based methods. 
Further discussion addresses the issue and possibilities for a moving grid based hydrocode invoking the DTFE, 
^ I and Delaunay tessellations, in an attempt to combine the virtues of the Eulerian and Lagrangian approaches. 

Key words, hydrodynamics, methods: N-body simulations, methods: numerical 

. >^ , 1. Introduction tinuous density field. The basic feature of the SPH pro- 
cedure for density estimation is based upon a convolution 
Smoothed Particle Hydrodynamics (SPH) has established tj^e discrete particle distribution with a particular user- 
itself as the workhorse for a variet y of astrophysical specified kernel function W. For a sample of N particles, 
fluid dyna mical computations (Lucy |1977| Ginghold & ^j^h masses ruj and locations r^, the density p at the lo- 
Monaghan[L977|). In a wide range of astrophysical envi- cation r, of particle i is given by 
ronments this Lagrangian scheme offers substantial and 
often crucial advantages over Eulerian, usually grid-based, 
schemes. Astrophysical applications such as cosmic struc- 
ture formation and galaxy formation, the dynamics of ■'^^ 

accretion disks and the formation of stars and plane- in which the kernel resolution is determined through the 

tary systems are testament to its versatility and sue- smoothing scale h,. Notice that generically the scale hi 

cesful performance (for an enumeration of applications, may be different for each individual particle, and thus may 

and corre spondi ng references, s ee e.g . the reviews by be set to adapt to the local particle density. Usually the 

MonaghanHnnaand BertschingerHnnHI)- functional dependence of the kernel W is chosen to be 

A crucial aspect of the SPH procedure concerns the spherically symmetric, so that it is a function of jr^ — 

proper estimation of the local density, i.e. the density at only. 

the location of the particles which are supposed to rep- The evolution of the physical system under considera- 

resent a fair - discrete - sampling of the underlying con- tion is fully determined by the movement of the discrete 
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particles. Given a properly defined density estimation pro- 
cedure, the equations of motion for the set of particles are 
specified through a suitable Lagrangian, if necessary in- 
cluding additional viscous forces (see e.g. Rasio [TMU|l . 

In practical implementations, however, the SPH pro- 
cedure involves a considerable number of artefacts. These 
stem from the fact that SPH particles represent functional 
averages over a certain Lagrangian volume. This averag- 
ing procedure is further aggravated by the fact that it 
is based upon a rather arbitrary user-specified choice of 
both the adopted resolution scale(s) hi and the functional 
form of the kernel W. Such a description of a physical 
system in terms of user-defined fuzzy clouds of matter is 
known to lead to considerable complications in realistic as- 
trophysical circumstances. Often, these environments in- 
volve fluid flows exhibiting complex spatial patterns and 
geometries. In particular in configurations characterized 
by strong gradients in physical characteristics - of which 
the density, pressure and temperature discontinuities in 
and around shock waves represent the most frequently en- 
countered example - SPH has been hindered by its relative 
inefficiency in resolving these gradients. 

Given the necessity for the user to specify the char- 
acteristics and parameter values of the density estimation 
procedure, the accuracy and adaptibility of the resulting 
SPH implementation hinges on the ability to resolve steep 
density contrasts and the capacity to adapt itself to the 
geometry and morphology of the local matter distribu- 
tion. A considerable improvement with respect to the early 
SPH implementations, which were based on a uniform 
smoothing length h, involves the use of adaptive smooth- 
ing lengths hi (Hernquist & Katz HQXl^ll . which provides 
the SPH calculations with a larger dynamic range and 
higher spatial resolution. The mass distribution in many 
(astro)physical systems and circumstances is often char- 
acterized by the presence of salient anisotropic patterns, 
usually identified as filamentary or planar features. To 
deal with such configurations, additional modifications in 
a few sophisticated implementations attempted to replace 
the conventional - and often unrealistic and restrictive - 
spherically symmetric kernels by ones whose configuration 
is more akin to the shape of the local mass distribution. 
The corresponding results do indeed represent a strong ar- 
gument for the importance of using geometrically adaptive 
density estimates. A noteworthy example is the introduc- 
tion of ellipsoidal kernels by Shapiro et al. Their 
shapes are stretched in accordance with the local flow. Yet, 
while evidently being conceptually superior, their practi- 
cal implementation does constitute a major obstacle and 
has prevented widescale use. This may be ascribed largely 
to the rapidly increasing number of degrees of freedom 
needed to specify and maintain the kernel properties dur- 
ing a simulation. 

Even despite their obvious benefits and improve- 
ments, these methods are all dependent upon the artificial 
parametrization of the local spatial density distribution in 
terms of the smoothing kernels. Moreover, the specifica- 
tion of the information on the density distribution in terms 



of extra non physical variables, necessary for the definition 
and evolution of the properties of the smoothing kernels, is 
often cumbersome to implement and may introduce subtle 
errors (Hernquist 110081 see however Nelson & Papaloizou 
ITHMI and Springel k Hernquist |2002 ). In many astrophys- 
ical applications this may lead to systematic artefacts in 
the outcome for the related physical phenomena. Within 
a cosmological context, for example, the X-ray visibility of 
clusters of galaxies is sensitively dependent upon the value 
of the local density, setting the intensity of the emitted X- 
ray emission by the hot intergalactic gas. This will be even 
more critical in the presence of feedback processes, which 
for sure will be playing a role when addressing the amount 
of predicted star formation in simulation studies of galaxy 
formation. 

Here, we seek to circumvent the complications induced 
by the kernel parametrization and introduce and propose 
an alternative to the use of kernels for the quantification of 
the density within the SPH formalism. This new method, 
based upon the Delaunay Tessellation Field Estimator 
(DTFE, Schaap & van de Weygaert 110001) , has been de- 
vised to mould and fully adapt itself to the configura- 
tion of the particle distribution. Unlike conventional SPH 
methods, it is able to deal self-consistently and naturally 
with anisotropics in the matter distribution, even when it 
concerns caustic-like transitions. In addition, it manages 
to succesfuUy treat density fields marked by structural 
features over a vast (dynamic) range of scales. 

The DTFE produces density estimates on the basis of 
the particle distribution, which is supposed to form a dis- 
crete spatial sampling of the underlying continuous den- 
sity field. As a linear multidimensional field interpolation 
algorithm it may be regarded as a first-order version of 
the natural neighbour algorithm for spatial interpolation 
(Sibson HOUl also see e.g. Okabe et al. I2000p . In gen- 
eral, applications of the DTFE to spatial point distribu- 
tions have demonstrated its success in dealing with the 
complications of anisotropic geometry and dynamic range 
(Schaap & van de Wevgaert l2l)Ul)|l . The key ingredient of 
the DTFE procedure is that of the Delaunay triangula- 
tion, serving as the complete covering of a sample volume 
by mutually exclusive multidimensional linear interpola- 
tion intervals. 

Delaunay tessellations (Delone 1934; see e.g. Okabe 
et al. 2000 for extensive review) form the natural frame- 
work in which to discuss the properties of discrete point 
sets, and thus also of discrete samplings of continuous 
fields. Their versatility and significance have been un- 
derlined by their widespread applications in such areas 
as computer graphics, geographical mapping and medi- 
cal imaging. Also, they have already found widespread 
application in a variety of 'conventional'grid-based fiuid 
dynamical computation schemes. This may concern their 
use as a non-regular application-oriented grid covering 
of physicalsystems, which represents a prominent proce- 
dure in technological applications. More innovating has 
been their use in Lagrangian 'moving-grid' implementa- 
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tions (see Mavripilis 119971 for a review, and Whitehurst 
EMI for a promising astrophysical application) . 

It seems therefore a good idea to explore the possibil- 
ities of applying the DTFE in the context of a numerical 
hydrodynamics code. Here, as a first step, we wish to ob- 
tain an idea of the performance of a hydro code involving 
the use of DTFE estimates with respect to an equivalent 
code involving regular SPH density estimates. The quality 
of the new DTFE method with respect to the conventional 
SPH estimates, and their advantages and disadvantages 
under various circumstances, are evaluated by a compari- 
son between the density field which would be yielded by a 
DTFE processing of the resulting SPH particle distribu- 
tion and that of the regular SPH procedure itself. In this 
study, we operate along these lines by a comparison of the 
resulting matter distributions in the situation of a repre- 
sentative stochastic multiphase density field. This allows 
us to make a comparison between both density estimates 
in a regime for which an improved method for density es- 
timates would be of great value. We should point out a 
major drawback of our approach, in that we do not re- 
ally treat the DTFE density estimate in a self-consistent 
fashion. Instead of being part of the dynamical equations 
themselves we only use it as an analysis tool of the pro- 
duced particle distribution. Nevertheless, it will still show 
the value of the DTFE in particle gasdynamics and give 
an indication of what kind of differences may be expected 
when incorporating in a fully self-consistent manner the 
DTFE estimate in an hydrocode. 

On the basis of our study, we will elaborate on the po- 
tential benefits of a hydrodynamics scheme based on the 
DTFE. Specifically, we outline how we would set out to 
develop a complete particle hydrodynamics code whose ar- 
tificial kernel based nature is replaced by the more natural 
and self-adaptive approach of the DTFE. Such a DTFE 
based particle hydrodynamics code would form a promis- 
ing step towards the development of a fully tessellation 
based quasi-Eulerian moving-grid hydrodynamical code. 
Such would yield a major and significant step towards 
defining a much needed alternative and complement to 
currently available simulation tools. 

2. DTFE and SPH density estimates 

The methods we use for SPH and DTFE density estimates 
have been extensively described elsewhere (Hernquist & 
KatzHnnni Schaap & van de Wevgaert IMIHI . Here, we 
will only summarize their main, and relevant, aspects. 

2.1. SPH density estimate 

Amongst the various density recepies employed within 
available SPH codes, we use the Hernquist & Katz l|1989|l 
symmetrized form of Eq. ^ using adaptive smoothing 
lengths: 

- ^^m, {W(|n-r,|,/i,) + VF(h-r,U,)}, (2) 
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The smoothing lengths hi are chosen such that the sum 
involves around 40 nearest neighbours. For the kernel 
W we take the conventional spline kernel described by 
Monaghan H1992|l . Other variants of the SPH estimate 
produce comparable results. 

2.2. DTFE density estimate 

The DTFE density estimating procedure consists of three 
basic steps. 

Starting from the sample of particle locations, the 
first step involves the computation of the correspond- 
ing Delaunay tessellation. Each Delaunay cell r,„ is the 
uniquely defined tetrahedron whose four vertices (in 3- 
D) are the set of 4 sample particles whose circumscribing 
sphere does not contain any of the other particles in the 
set. The Delaunay tessellation is the full covering of space 
by the complete set of these mutually disjunct tetrahedra. 
Delaunay tessellations are well known concepts in stochas- 
tic and computational geometry (Delaunay .1934. for fur- 
ther references see e.g. Okabe et al. l2()()()l M(Ziller [r994l and 
van de Wevgaert ll991|l . 

The second step involves estimating the density at the 
location of each of the particles in the sample. From the 
definition of the Delaunay tessellation, it may be evident 
that there is a close relationship between the volume of a 
Delaunay tetrahedron and the local density of the gen- 
erating point process (telling examples of this may be 
seen in e.g. Schaap & van de Weygaert 2002a). Evidently, 
the "empty" cirumscribing spheres corresponding to the 
Delaunay tetrahedra, and the volumes of the resulting 
Delaunay tetrahedra, will be smaller as the number den- 
sity of sample points increases, and vice versa. Following 
this observation, a proper density estimate p at the loca- 
tion Xi of a sampling point i is obtained by determining 
the properly calibrated inverse of the volume Wvor.i of 
the corresponding contiguous Voronoi cell. The contigu- 
ous Voronoi cell Wvor,i is the union of all Delaunay tetra- 
hedra Tm.i of which the particle i forms one of the four 
vertices, i.e. Wvor,i = Um^m,i- general, when a parti- 
cle i is surrounded by Nt Delaunay tetrahedra, each with 
a volume V{Tm^i), the volume of the resulting contiguous 
Voronoi cell is 

Nt 
'm=l 

Note that Nt is not a constant, but in general may acquire 
a different value for each point in the sample. For a Poisson 
distribution of particles this is a non-integer number in the 
order of (Nt) ~ 27 (van de Wevgaert ll994|l . Generalizing 
to an arbitrary I?-dimensional space, and assuming that 
each particle i has been assigned a mass m^, the estimated 
density pi at the location of particle i is given by (see 
Schaap & van de Weygaert I2000|l 

pir,) = (D + l)-^, (4) 
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SPH DTFE 

Fig. 1. Comparison of the DTFE performance versus that of the regular SPH method in a characteristic configuration, 
that of a hydrodynamic simulation of the multiphase interstellar medium. Top left panel: the particle distribution in 
a 0.6 X 0.6 kpc simulation region, within a slice with a width of 0.005 kpc. Bottom left frame: 2-D slice through the 
resulting (3-D) SPH density field reconstruction. Bottom right frame: the corresponding (3-D) density field recon- 
struction produced by the DTFE procedure. Top righthand frame: summary, in terms of a quantitative point-by-point 
comparison between the DTFE and SPH density estimates, pdtfe and psPH- Abscissa: the value of the SPH density 
estimate (normalized by the average density (p)). Ordinate: the ratio of DTFE estimate to the SPH density estimate, 
Pdtfe/psph- These quantities are plotted for each particle location in the full simulation box. 



In this, we explicitly express Wvor,i for the general D- 
dimensional case. The factor {D + 1) is a normaliza- 
tion factor, accounting for the {D + 1) different con- 
tiguous Voronoi hypercells to which each Delaunay hy- 
per "tetrahedron" is assigned, one for each vertex of a 
Delaunay hyper "tetrahedron" . 



The third step is the interpolation of the estimated 
densities pi over the full sample volume. In this, the DTFE 
bases itself upon the fact that each Delaunay tetrahedron 
may be considered the natural multidimensional equiva- 
lent of a linear interpolation interval (see e.g Bernardeau 
& van de Weveaert 11996(1 . Given the {D + 1) vertices of 
a Delaunay tetrahedron with corresponding density es- 
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timates pj, the value p(r) at any location r within the 
tetrahedron can be straightforwardly determined by sim- 
ple linear interpolation, 

p{y) = /5(rjo) + (Vp)DoKm • (r - r^o) , (5) 

in which r,;o is the location of one of the Delaunay vertices 
i. This is a trivial evaluation once the value of the (linear) 
density gradient (Vp)Dci.m has been estimated. For each 
Delaunay tetrahedron this is accomplished by solv- 
ing the the system of D linear equations corresponding to 
each of the remaining D Delaunay vertices constituting 
the Delaunay tetrahedron Tm- The "minimum triangula- 
tion" property of Delaunay tessellations underlying this 
linear interpolation, minimum in the sense of representing 
a volume-covering network of optimally compact multi- 
dimensional "triangles" , has been a well-known property 
utilized in a variety of imaging and surface rendering ap- 
plications such as geographical mapping and various com- 
puter imaging algorithms. 

2.3. Comparison 

Comparing the two methods, we see that in the case of 
SPH the particle 'size' and 'shape' (i.e. its domain of in- 
fluence) is determined by some arbitrary kernel W{r^ hi) 
and a fortuitous choice of smoothing length hi (assuming, 
along with the major share of SPH procedures, a radially 
symmetric kernel). In the case of the DTFE method the 
particles' influence region is fully determined by the sizes 
and shapes of the Delaunay cells Tm,i, themselves solely 
dependent on the particle distribution. In other words, in 
regular SPH the density is determined through the ker- 
nel function W{x), while in DTFE it is solely the parti- 
cle distribution itself setting the estimated values of the 
density. Contrary to the generic situation for the kernel 
dependent methods, there are no extra variables left to be 
determined. One major additional advantage is that it is 
therefore not necessary to worry about the evolution of 
the kernel parameters. 

Both methods do display some characteristic artefacts 
in their density reconstructions (see Fig. 1). To a large ex- 
tent these may be traced back to the implicit assumptions 
involved in the interpolation procedures, a necessary con- 
sequence of the finite amount of information contained in 
a discrete representation of a continuous field. SPH den- 
sity fields implicitly contain the imprint of the specified 
and applied kernel which, as has been discussed before, 
may seriously impart its resolving power and capacity to 
trace the true geometry of structures. The DTFE tech- 
nique, on the other hand, does produce triangular arte- 
facts. At instances conspicuously visible in the DTFE re- 
constructed density fields, they are the result of the linear 
interpolation scheme employed for the density estimation 
at the locations not coinciding with the particle positions. 
In principle, this may be substantially improved by the 
use of higher order interpolation schemes. Such higher- 
order schemes have indeed been developed, and the ones 



based upon the natural neighbour interpolation prescrip- 
tion of Sibson (,1981) have already been succesfuUy applied 
to two-dimensional problems in the field of geophysics 
(Sambridge et al. 119951 Braun & Sambridge I1995|l and 
solid state physics fSukiimar ll998|l . 

3. Case study: two-phase interstellar medium 

For the sake of testing and comparing the SPH and DTFE 
methods, we assess a snapshot from a simulation of the 
neutral ISM. The model of the ISM is chosen as an illus- 
tration rather than as a realistic model. 

The "simulation" sample of the ISM consists of HI 
gas confined in a periodic simulation box with a size 
L = 0.6 kpc"^. The initially uniform density of the gas 
is jiH ~ 0.3 cm~^, while its temperature is taken to be 
T — 10000 K. No fluctuation spectrum is imposed to set 
the initial featureless spatial gas distribution. To set the 
corresponding initial spatial distribution of the N = 64000 
simulation particles, we start from relaxed initial condi- 
tions according to a "glass" distribution (e.g. White IT^MIl . 

The evolution of the gas is solely a consequence of fluid 
dynamical and thcrmodynamical processes. No self gravity 
is included. As for the thcrmodynamical state of the gas, 
cooling is implemented using a fit to the Dalgarno-McCray 
l|1972|l cooling curve. The heating of the gas is accom- 
plished through photo-electric grain heating, attributed to 
a constant FUV background (1.7 Go, with Go the Habing 
field) radiation field. The parameters are chosen such that 
after about 15 Myrs a two-phase medium forms which 
consists of warm (10000 K) and cold (> 100 K) HI gas. 

The stage at which a two-phase medium emerges forms 
a suitable point to investigate the performance of the SPH 
and DTFE methods. At this stage we took a snapshot 
from the simulation, and subjected it to further analysis. 
For a variety of reasons, the spatial gas distribution of the 
snapshot is expected to represent a challenging configu- 
ration. The multiphase character of the resulting particle 
configuration is likely to present a problem for regular 
SPH. Density contrasts of about four orders of magni- 
tude separate dense clumps from the surrounding diffuse 
medium through which they are dispersed. Note that a 
failure to recover the correct density may have serious 
repercussions for the computed effects of cooling. In ad- 
dition, we notice the presence of physical structures with 
conspicuous, aspherical geometries (see Fig. 1 & 2), such 
as anisotropic sheets and filaments as well as dense and 
compact clumps, which certainly do form a challenging 
aspect for the different methods. 

3.1. Results 

Fig- n offers a visual impression of the differences in per- 
formance between the SPH and DTFE density reconstruc- 
tions. The greyscale density maps in Fig. 1 (lower left: 
SPH, lower right: DTFE) represent 2-D cuts through the 
corresponding 3-D density field reconstructions (note that 
contrary to the finite width of the corresponding particle 
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Fig. 2. Systematic analysis of the differences between the DTFE and SPH density estimates, pdtfe and psPH- Basis 
of the analysis is a point-by-point comparison of these two density estimates. Top lefthand frame: diagram of the 
value of the ratio pdtfe/psph (ordinate) versus psph/(p) (abscissa) for each of the points in the simulation volume. 
Indicated in this scatter diagram are four sectors, each of which corresponds to particles residing in a physically 
different regime/phase. On the basis of this identification, the full set of particles is dissected into the corresponding 
four composing particle samples. Top righthand frame: the spatial distribution of the full set of particles in a 0.04 kpc 
wide slice. The subsequent 4 frames (from central left to bottom right) show, for each indicated sector in the scatter 
diagram, the spatial distribution of the corresponding particles (within the same 0.04 kpc slice). 



slice, upper left frame, these constitute planes with zero 
thickness). 

Immediately visible is the more crispy appearance of 
the DTFE density field, displaying substantially more con- 
trast in conjunction with more pronounced structural fea- 
tures. Look e.g. at the compact clump in the lower right- 
hand corner {X » 0.5, y « 0.12), forming a prominent 
and tight spot in the DTFE density field. The clump at 
{X « 0.48, Y « 0.52) represents another telling example, 
visible as a striking peak in the DTFE rendering while 
hardly noticeable in the SPH reconstruction. Structures 
in the SPH field have a more extended appearance than 
their counterparts in the DTFE field, whose matter con- 
tent has been smeared out more evenly, over a larger vol- 
ume, yielding features with a significantly lower contrast. 
In this assessment it becomes clear that the DTFE re- 
construction adheres considerably closer to the original 
particle distribution (top lefthand frame). Apparently the 
DTFE succeeds better in rendering the shapes, the coher- 
ence and the internal composition in the displayed particle 
distribution. At various locations, the DTFE even man- 
ages to capture structural details which seem to be absent 
in the SPH density field. 

To quantify the visual impressions of Fig. 1, and to 
analyze the nature of the differences between the two 
methods, we plot the ratio pdtfe/psph as a function of 
the SPH density estimate Psph/(p) (in units of the aver- 
age density (p)). Doing so for all particles in the sample 
(Fig. 1, top righthand. Fig. 2, top lefthand) immediately 
reveals interesting behaviour. The scatter diagram does 
show that the discrepancies between the two methods may 
be substantial, with density estimates at various instances 
differing by a factor of 5 or more. 

Most interesting is the finding that we may distinguish 
clearly identifiable and distinct regimes in the scatter dia- 
gram of pdtfe/psph versus pspa/ip)- Four different sec- 
tors may be identified in the scatter diagram. Allowing for 
some arbitrariness in their definition, and indicating these 
regions by digits 1 to 4, we may organize the particles ac- 
cording to density-related criteria, roughly specified as (we 
refer to Fig. 2, top left frame, for the precise definitions of 
the domains): 

1. low density regions: 

Psph/(p) < 1 

2. medium density regions, DTFE smaller than SPH: 

Pdtfe < Psph; 1 < Pspn/ip) < 10 



3. medium density regions, DTFE larger than SPH: 

Pdtfe > Psph; 1 < Psph/(p) < 10 

4. high density regions: 

Pdtfe > Psph; Psph/(p) > 10 

The physical meaning of the distinct sectors in the scat- 
ter diagram becomes apparent when relating the various 
regimes with the spatial distribution of the corresponding 
particles. This may be appreciated from the five subse- 
quent frames in Fig. 2, each depicting the related particle 
distribution in the same slice of width 0.04 kpc. The cen- 
tre and bottom frames, numbered 1 to 4, show the spatial 
distribution of each group of particles, isolated from the 
complete distribution (top right frame, Fig. 2). These par- 
ticle slices immediately reveal the close correspondence be- 
tween any of the sectors in the scatter diagram and typical 
features in the spatial matter distribution of the two-phase 
interstellar medium. This systematic behaviour seems to 
point to truly fundamental differences in the workings of 
the SPH and DTFE methods, and would be hard to un- 
derstand in terms of random errors. The separate spatial 
features in the gas distribution seem to react differently 
to the use of the DTFE method. 

We argue that the major share of the disparity be- 
tween the SPH and DTFE density estimates has to be 
attributed to SPH, mainly on the grounds of the known 
fact that SPH is poor in handling nontrivial configura- 
tions such as encountered in multiphase media. By sepa- 
rately assessing each regime, we may come to appreciate 
how these differences arise. In sector 1, involving the dif- 
fuse low density medium, the DTFE and SPH estimates 
are of comparable magnitude, be it that we do observe a 
systematic tendency. In the lowest density realms, whose 
relatively smooth density does not raise serious obstacles 
for either method, DTFE and SPH are indeed equal (with 
the exception of variations to be attributed to random 
noise). However, near the edges of the low density regions, 
SPH starts to overestimate the local density as the ker- 
nels do include particles within the surrounding high den- 
sity structures. The geometric interpolation of the DTFE 
manages to avoid this systematic effect (see e.g. Schaap & 
van de Weygaert 2002a, b), which explains the systematic 
linear decrease of the ratio Pdtfe/psph with increasing 
Psph/(p)- To the other extreme, the high density regions 
in sector 4 are identified with compact dense clumps as 
well as with their extensions into connecting filaments 
and walls. On average DTFE yields higher density esti- 
mates than SPH, frequently displaying superior spatial 
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resolution (see also greyscale plot in Fig. 1). Note that 
the repercussions may be far-reaching in the context of a 
wide variety of astrophysical environments characterized 
by strongly density dependent physical phenomena and 
processes ! The intermediate regime of sectors 2 and 3 
clearly connects to the filamentary structures in the gas 
distribution. Sector 2, in which the DTFE estimates are 
larger than those of SPH, appears to select out the inner 
parts of the filaments and walls. By contrast, the higher 
values for the SPH produced densities in sector 3 are re- 
lated to the outer realms of these features. This charac- 
teristic distinction can be traced back to the failure of 
the SPH procedure to cope with highly anisotropic parti- 
cle configurations. While it attempts to maintain a fixed 
number of neighbours within a spherical kernel, it smears 
out the density in a direction perpendicular to the fila- 
ment. This produces lower estimates in the central parts, 
which are compensated for with higher estimates in the 
periphery. Evidently, the adaptive nature of DTFE does 
not appear to produce similar deficiencies. 

4. The DTFE particle method 

Having demonstrated the improvement in quality of the 
DTFE density estimates, this suggests a considerable po- 
tential for incorporating the DTFE in a self-consistent 
manner within a hydrodynamical code. Here, we first wish 
to indicate a possible route for accomplishing this in a 
particle hydrodynamics code through replacement of the 
kernel based density estimates by the DTFE density 
estimates. We are currently in the process of implement- 
ing this. The formalism on which this implementation is 
based can be easily derived, involving nontrivial yet mi- 
nor modifications. Essentially, it uses the same dynamic 
equations for gas particles as those in the regular SPH for- 
malism, the fundamental adjustment being the insertion 
of the DTFE densities instead of the regular SPH ones. In 
addition, a further difference may be introduced through 
a change in treatment of viscous forces. Ultimately, this 
will work out into different equations of motion for the 
gas particles. A fundamental property of a DTFE based 
hydrocode, by construction, is that it conserves mass ex- 
actly and therefore obeys the continuity equation. This is 
not necessarily true for SPH implementations (Hernquist 

& KatzCnHni)- 

The start of the suggested DTFE particle method is 
formed by the discretized expression for the Lagrangian L 
for a compressible, nondissipative flow. 



(6) 



where rrii is the mass of particle i, Vi its velocity, Si the 
corresponding entropy and m its specific internal energy. 
In this expression, pi is the density at location z, as yet 
unspecified. The resulting Euler-Lagrange equations are 



dt 



3 



dxi 



(7) 



The standard SPH equations of motion then follow after 
inserting the SPH density estimate (Eq.^. Instead, inser- 
tion of the DTFE density (Eq. IH will lead to the corre- 
sponding equations of motion for the DTFE-based formal- 
ism. Note that the usual conservation properties related to 
Eq. 1^ remain intact. After some algebraic manipulation, 
thereby using the basic thermodynamic relation for a gas 
with equation of state 



dpi 



El 

^2 



(8) 



we finally obtain the equations of motion for the gas par- 
ticles (moving in D-dimensional space), 



dv 1 

" 'd + i ^ ^^^"''^ 



dxi 



(9) 



This expression involves a summation over all Nt 
Delaunay tetrahedra TVi.i, with volumes V{Tm^^)^ which 
have the particle i as one of its four vertices. The pressure 
term P{Tm,i) is the sum over the pressures Pj at the four 
vertices j of tetrahedron Tm,i, P{Tm,i) = Y^Pj- 

As an interesting aside, we point out that unlike in 
the conventional SPH formalism, this procedure implies 
an exactly vanishing acceleration dvi/dt in the case of a 
constant pressure P at each of the vertices of the Delaunay 
tetrahedra containing particle i as one of their vertices. 
The reason for this is that one can then invoke the defi- 
nition of the volume of the contiguous Voronoi cell corre- 
sponding to point i (Eqn. yielding 



dvt 
'dt 



1 



D + 1 



P 



dWv 



dxi 



(10) 



Since the volume of the contiguous Voronoi cell does not 
depend on the position of particle i itself (it lies in the 
interior of the contiguous Voronoi cell), the resulting ac- 
celeration vanishes. Another interesting notion, which was 
pointed out by Icke H2U(J2|I . is that Delaunay tessellations 
also provide a unique opportunity to include a natural 
treatment of the viscous stresses in the physical system. 
We intend to elaborate on this possibility in subsequent 
work dealing with the practical implementation along the 
lines sketched above. 

5. Delaunay tessellations 

and 'moving grid' hydrocodes 

Ultimately, the ideal hydrodynamical code would combine 
the advantages of the Eulerian as well as of the Lagrangian 
approach. In their simplest formulation, Eulerian algo- 
rithms cover the volume of study with a fixed grid and 
compute the fluid transfer through the faces of the (fixed) 
grid cell volumes to follow the evolution of the system. 
Lagrangian formulations, on the other hand, compute the 
system by following the ever changing volume and shape 
of a particular individual element of gas (interestingly, the 
'Lagrangian' formulation is also due to Euler 118621 who 
employed this formalism in a letter to Lagrange, who later 
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proposed these ideas in a publication by himself. [17621 see 
Whitehurstll995). 

For a substantial part the success of the DTFE may 
be ascribed to the use of Delaunay tessellations as an 
optimally covering grid. This suggests that they may 
also be ideal for the use in moving grid implementations 
for hydrodynamical calculations. As in our SPH appli- 
cation, such hydrocodes with Delaunay tessellations at 
their core would warrant a close connection to the un- 
derlying matter distribution. Indeed, attempts towards 
such implementations have already been introduced in 
the context of a few specific, mainly two-dimensional, ap- 
plications (Whitehurst 119951 Braun & Sambridge .1995. 
Sukumar 11998(1 . Alternative attempts towards the devel- 
opment of moving grid codes, in an astrophysical context, 
have shown their potential (Gnedin il995i Pen iil998p . 

For a variety of astrophysical problems it is indeed es- 
sential to have such advanced codes at one's disposal. An 
example of high current interest may offer a good illustra- 
tion. Such an example is the reionization of the intergalac- 
tic medium by the ionizing radiation emitted by the first 
generation of stars, (proto)galaxies and/or active galactic 
nuclei. These radiation sources will form in the densest 
regions of the universe. To be able to resolve these in suf- 
ficient detail, it is crucial that the code is able to focus 
in onto these densest spots. Their emphasis on mass res- 
olution makes Lagrangian codes - including SPH - usu- 
ally better equipped to do so, be it not yet optimally. 
On the other hand, it is in the low density regions that 
most radiation is absorbed at first. In the early stages the 
reionization process is therefore restricted to the huge un- 
derdense fraction of space. Simulation codes should there- 
fore properly represent and resolve the gas density distri- 
bution within these voidlike regions. The uniform spatial 
resolution of the Eulcrian codes is better suited to ac- 
complish this. Ideally, however, a simulation code should 
be able to combine the virtues of both approaches, yield- 
ing optimal mass resolution in the high density source re- 
gions and a proper coverage of the large underdense re- 
gions. Moving grid methods, of which Delaunay tessella- 
tion based ones will be a natural example, may indeed 
be the best alternative, as the reionization simulations by 
Gnedin l,1995|i) appear to indicate. There have been many 
efforts in the context of Eulerian codes towards the devel- 
opment of Adaptive Mesh Refinement ( AMR) algorithms ( 
Berger .1989)) , which have achieved a degree of maturity. 
Their chief advantage is their ability to concentrate com- 
putational effort on regions based on arbitrary refinement 
criteria, where, in the basic form at least, moving grid 
methods refine on a mass resolution criterion. However 
they are still constrained by the use of regular grids, which 
may introduce artifacts due to the presence of preferred 
directions in the grid. The advantages of a moving grid 
fluid dynamics code based on Delaunay tessellations have 
been most explicitly demonstrated by the implementation 
of a two-dimensional lagrangian hydrocode (FLAME) by 
Whitehurst H1995|l . These advantages will in principle ap- 
ply to any such algorithm, in particular also for three- 



dimensional implementations (of which we are currently 
unaware). Whitehurst H1995|l enumerated various poten- 
tial benefits in comparison with conventional SPH codes, 
most importantly the following: 

1. SPH needs a smoothing length h. 

2. SPH needs an arbitrary kernel function W. 

3. The moving grid method does not need an (unphysical) 
artificial viscosity to stabilize solutions. 

The validity of the first two claims has of course also been 
demonstrated in this study for particle methods based on 
DTFE. Whitehurst showed additionally that there is an 
advantage of moving grid methods over Eulerian grid- 
based ones. The implementation of Whitehurst, which 
used a first-order solver and a limit on the shape of grid 
cells to control the effects of shearing of the grid, was far 
superior to all tested first-order Eulerian codes, and su- 
perior to many second-order ones as well. The adaptive 
nature of the Langrangian method and the fact that the 
resulting grid has no preferred directions are key factors in 
determining the performance of moving grid methods such 
as FLAME. For additional convincing arguments, includ- 
ing the other claims, we may refer the reader to the truly 
impressive case studies presented by Whitehurst H1995|l . 

6. Summary &. discussion 

Here we have introduced the DTFE as an alternative den- 
sity estimator for particle fluid dynamics. Its principle as- 
set is that it is fully self-adaptive, resulting in a density 
field reconstruction which closely reproduces, usually in 
meticulous detail, the characteristics of the spatial particle 
distribution. It may do so because of its complete indepen- 
dence of arbitrary user-specified smoothing functions and 
parameters. Unlike conventional methods, such as the ker- 
nel estimators used in SPH, it manages to faithfully repro- 
duce the anisotropics in the local particle distribution. It 
therefore automatically refiects the genuine geometry and 
shape of the structures present in the underlying density 
field. This is in marked contrast with kernel based meth- 
ods, which almost without exception produce distorted 
shapes of density features, the result of the convolution of 
the real structure with the intrinsic shape of the smooth- 
ing function. Its adaptive and local nature also makes it 
optimally suited for reconstructing the hierarchy of scales 
present in the density distribution. In kernel based meth- 
ods the internal structural richness of density features is 
usually suppressed on scales below that of the character- 
istic (local) kernel scale. DTFE, however, is solely based 
upon the particle distribution itself and follows the density 
field wherever the discrete representation by the particle 
distribution allows it to do so. Its capacity to resolve struc- 
tures over a large dynamic range may prove to be highly 
beneficial in many astrophysical circumstances, quite of- 
ten involving environments in which we encounter a hier- 
archical embedding of small-scale structures within more 
extended ones. 
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In this study we have investigated the performance of 
the DTFE density estimator in the context of a Smooth 
Particle Hydrodynamics simulation of a multiphase in- 
terstellar medium of neutral gas. The limited spatial 
resolution of current particle hydrodynamics codes are 
known to implicate considerable problems near regions 
with e.g. steep density and temperature gradients. In par- 
ticular their handling of shocks forms a source of consider- 
able concern. SPH often fails in and around these regions, 
so often playing a critical and vital role in the evolution of 
a physical system. Our study consists of a comparison and 
confrontation of the conventional SPH kernel based den- 
sity estimation procedure with the corresponding DTFE 
density field reconstruction method. 

The comparison of the density field reconstructions 
demonstrated convincingly the considerable improvement 
embodied by the DTFE procedure. This is in particular 
true at locations and under conditions where SPH ap- 
pears to fail. Filamentary and sheetlike structures pro- 
vide telling examples of the superior DTFE handling with 
respect to the regular SPH method, with the most pro- 
nounced improvement occurring in the direction of the 
steepest density gradient. 

Having shown the success of the DTFE, we are con- 
vinced that its application towards the analysis of the out- 
come of SPH simulations will prove to be highly beneficial. 
This may be underlined by considering a fitting illustra- 
tion. Simulations of the settling and evolution of the X-ray 
emitting hot intracluster gas in forming clusters of galax- 
ies do represent an important and cosmologically relevant 
example (see Borgani & Guzzo El)(Jll and Rosati et al. 120021 
for recent reviews). The X-ray luminosity is strongly de- 
pendent upon the density of the gas. The poor accuracy 
of the density determination in regular SPH calculations 
therefore yields deficient X-ray luminosity estimates (see 
Bertschinger 119981 and Rosati et al. 120021 for relevant re- 
cent reviews). Despite a number of suggested remedies, 
such as separating particles according to their tempera- 
ture, their ad-hoc nature does not evoke a strong sense 
of confidence in the results. Numerical limitations will of 
course always imply a degree of artificial smoothing, but 
by invoking tools based upon the DTFE technique there is 
at least a guarantee of an optimal retrieval of information 
contained in the data. 

Despite its promise for the use in a variety of analysis 
tools for discrete data samples, such as particle distribu- 
tions in computer simulations or galaxy catalogues in an 
observational context, its potential would be most opti- 
mally exploited by building it into genuine new fluid dy- 
namics codes. Some specific (two-dimensional) examples 
of succesful attempts in other scientific fields were men- 
tioned, and we argue for a similar strategy in astrophysics. 
One path may be the upgrade of current particle hydrody- 
namics codes by inserting DTFE technology. In this study, 
we have outlined the development of such a SPH-like hy- 
drodynamics scheme in which the regular kernel estimates 
are replaced by DTFE estimates. One could interpret this 
in terms of the replacement of the user-specified kernel 



by the self-adaptive contiguous Delaunay cell, solely de- 
pendent on the local particle configuration. An additional 
benefit will be that on the basis of the localized connec- 
tions in a Delaunay tessellations it will be possible to de- 
fine a more physically motivated artificial viscosity term. 

The ultimate hydrodynamics algorithm would com- 
bine the virtues of Eulerian and Lagrangian techniques. 
Considering the positive experiences with DTFE, it ap- 
pears to be worthwhile within the context of 'moving 
grid' or 'Lagrangian grid' methods to investigate the use 
of Delaunay tessellations for solving the Euler equations. 
With respect to a particle hydrodynamics code, the self- 
adaptive virtues of DTFE and its ability to handle ar- 
bitrary density jumps with only one intermediate point 
may lead to significant improvements in the resolution 
and shock handling properties. Yet, for grid based meth- 
ods major complications may be expected in dealing with 
the non-regular nature of the corresponding cells, compli- 
cating the handling of flux transport along the boundaries 
of the Delaunay tetrahedra. 

The computational cost of DTFE resembling tech- 
niques is not overriding. The CPU time necessary for gen- 
erating the Delaunay tessellation corresponding to a point 
set of N particles is in the order of O(A^logA^), compa- 
rable to the cost of generating the neighbour list in SPH. 
Within an evolving point distribution these tessellation 
construction procedures may be made far more efficient, 
as small steps in the development in the system will induce 
a correspondingly small number of tetrahedron (identity) 
changes. Such dynamic upgrading routines are presently 
under development. 

In summary, in this work we have argued for and 
demonstrated the potential and promise of a natural com- 
putational technique which is based upon one of the most 
fundamental and natural tilings of space, the Delaunay 
tessellation. Although the practical implementation will 
undoubtedly encounter a variety of complications, depen- 
dent upon the physical setting and scope of the code, 
the final benefit of a natural moving grid hydrodynamics 
code for a large number of astrophysical issues may not 
only represent a large progress in a computational sense. 
Its major significance may be found in its ability to ad- 
dress fundamental astrophysical problems in a new and 
truely natural way, leading to important new insights in 
the workings of the cosmos. 
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